﻿using System;
using System.Globalization;
using System.Text.RegularExpressions;

namespace QAToolsAddin {
    //http://home.hiwaay.net/~taylorc/toolbox/geography/geoutm.html
    //https://www.ibm.com/developerworks/java/library/j-coordconvert/

    class Coordinates {
        //http://www.jstott.me.uk/jscoord/

        public class LL {
            public double Latitude { get; private set; }
            public double Longitude { get; private set; }

            public static bool TryParseLL(object text, out object result) {
                result = null;
                if (text is string == false || string.IsNullOrEmpty((string)text)) return false;
                var split = ((string)text).Trim().Split(' ');
                double xx, yy;
                if (!(split.Length == 2
                    && double.TryParse(split[0].Trim(), NumberStyles.Float, CultureInfo.InvariantCulture, out xx)
                    && double.TryParse(split[1].Trim(), NumberStyles.Float, CultureInfo.InvariantCulture, out yy)))
                    return false;
                try {
                    result = new LL(xx, yy);
                    return true;
                } catch { }
                return false;
            }

            public static bool TryParseDMS(object text, out object result) {
                result = null;
                DMS dms;
                bool succeed = DMS.TryParse(text, out dms);
                if (succeed)
                    result = dms.ToLL();
                return succeed;
            }

            public static bool TryParseUTM(object text, out object result) {
                result = null;
                UTM utm;
                bool succeed = UTM.TryParse((string)text, out utm);
                if (succeed)
                    result = utm.ToLL();
                return succeed;
            }

            public LL(double latitude, double longitude) {
                this.Latitude = latitude;
                this.Longitude = longitude;
            }

            public double Distance(LL to) {
                double er = 6366.707;
                double latFrom = Deg2rad(this.Latitude);
                double latTo = Deg2rad(to.Latitude);
                double lngFrom = Deg2rad(this.Longitude);
                double lngTo = Deg2rad(to.Longitude);
                double x1 = er * Math.Cos(lngFrom) * Math.Sin(latFrom);
                double y1 = er * Math.Sin(lngFrom) * Math.Sin(latFrom);
                double z1 = er * Math.Cos(latFrom);
                double x2 = er * Math.Cos(lngTo) * Math.Sin(latTo);
                double y2 = er * Math.Sin(lngTo) * Math.Sin(latTo);
                double z2 = er * Math.Cos(latTo);
                return Math.Sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2) + (z1 - z2) * (z1 - z2));
            }

            public OSRef ToOSRef() {
                var airy1830 = new RefEll(6377563.396, 6356256.909);
                double OSGB_F0 = 0.9996012717;
                double N0 = -100000.0;
                double E0 = 400000.0;
                double phi0 = Deg2rad(49.0);
                double lambda0 = Deg2rad(-2.0);
                double a = airy1830.maj;
                double b = airy1830.min;
                double eSquared = airy1830.ecc;
                double phi = Deg2rad(this.Latitude);
                double lambda = Deg2rad(this.Longitude);
                double n = (a - b) / (a + b);
                double v = a * OSGB_F0 * Math.Pow(1.0 - eSquared * SinSquared(phi), -0.5);
                double rho = a * OSGB_F0 * (1.0 - eSquared) * Math.Pow(1.0 - eSquared * SinSquared(phi), -1.5);
                double etaSquared = (v / rho) - 1.0;
                double M = (b * OSGB_F0)
                    * (((1 + n + ((5.0 / 4.0) * n * n) + ((5.0 / 4.0) * n * n * n))
                    * (phi - phi0))
                    - (((3 * n) + (3 * n * n) + ((21.0 / 8.0) * n * n * n)) * Math.Sin(phi - phi0) * Math.Cos(phi + phi0))
                    + ((((15.0 / 8.0) * n * n) + ((15.0 / 8.0) * n * n * n)) * Math.Sin(2.0 * (phi - phi0)) * Math.Cos(2.0 * (phi + phi0)))
                    - (((35.0 / 24.0) * n * n * n) * Math.Sin(3.0 * (phi - phi0)) * Math.Cos(3.0 * (phi + phi0))));
                double I = M + N0;
                double II = (v / 2.0) * Math.Sin(phi) * Math.Cos(phi);
                double III = (v / 24.0)
                    * Math.Sin(phi)
                    * Math.Pow(Math.Cos(phi), 3.0)
                    * (5.0 - TanSquared(phi) + (9.0 * etaSquared));
                double IIIA = (v / 720.0)
                    * Math.Sin(phi)
                    * Math.Pow(Math.Cos(phi), 5.0)
                    * (61.0 - (58.0 * TanSquared(phi)) + Math.Pow(Math.Tan(phi), 4.0));
                double IV = v * Math.Cos(phi);
                double V = (v / 6.0) * Math.Pow(Math.Cos(phi), 3.0) * ((v / rho) - TanSquared(phi));
                double VI = (v / 120.0)
                    * Math.Pow(Math.Cos(phi), 5.0)
                    * (5.0 - (18.0 * TanSquared(phi)) + (Math.Pow(Math.Tan(phi), 4.0)) + (14 * etaSquared) - (58 * TanSquared(phi) * etaSquared));

                double N = I + (II * Math.Pow(lambda - lambda0, 2.0)) + (III * Math.Pow(lambda - lambda0, 4.0)) + (IIIA * Math.Pow(lambda - lambda0, 6.0));
                double E = E0 + (IV * (lambda - lambda0)) + (V * Math.Pow(lambda - lambda0, 3.0)) + (VI * Math.Pow(lambda - lambda0, 5.0));

                return new OSRef(E, N);
            }

            public UTM ToUTM() {
                var wgs84 = new RefEll(6378137, 6356752.314);
                double UTM_F0 = 0.9996;
                double a = wgs84.maj;
                double eSquared = wgs84.ecc;
                double longitude = this.Longitude;
                double latitude = this.Latitude;
                double latitudeRad = latitude * (Math.PI / 180.0);
                double longitudeRad = longitude * (Math.PI / 180.0);
                int longitudeZone = (int)Math.Floor((longitude + 180.0) / 6.0) + 1;

                // Special zone for Norway
                if (latitude >= 56.0 && latitude < 64.0 && longitude >= 3.0 && longitude < 12.0)
                    longitudeZone = 32;

                // Special zones for Svalbard
                if (latitude >= 72.0 && latitude < 84.0) {
                    if (longitude >= 0.0 && longitude < 9.0) longitudeZone = 31;
                    else if (longitude >= 9.0 && longitude < 21.0) longitudeZone = 33;
                    else if (longitude >= 21.0 && longitude < 33.0) longitudeZone = 35;
                    else if (longitude >= 33.0 && longitude < 42.0) longitudeZone = 37;
                }

                double longitudeOrigin = (longitudeZone - 1.0) * 6.0 - 180.0 + 3.0;
                double longitudeOriginRad = longitudeOrigin * (Math.PI / 180.0);
                char UTMZone = GetUTMLatitudeZoneLetter(latitude);
                double ePrimeSquared = (eSquared) / (1.0 - eSquared);
                double n = a / Math.Sqrt(1 - eSquared * Math.Sin(latitudeRad) * Math.Sin(latitudeRad));
                double t = Math.Tan(latitudeRad) * Math.Tan(latitudeRad);
                double c = ePrimeSquared * Math.Cos(latitudeRad) * Math.Cos(latitudeRad);
                double A = Math.Cos(latitudeRad) * (longitudeRad - longitudeOriginRad);
                double eSquared2 = eSquared * eSquared;
                double eSquared3 = eSquared * eSquared * eSquared;

                double M = a * (
                    (1.0 - eSquared / 4.0 - 3.0 * eSquared2 / 64.0 - 5.0 * eSquared3 / 256.0) * latitudeRad
                    - (3.0 * eSquared / 8.0 + 3.0 * eSquared2 / 32.0 + 45.0 * eSquared3 / 1024.0) * Math.Sin(2.0 * latitudeRad)
                    + (15.0 * eSquared * eSquared / 256.0 + 45.0 * eSquared3 / 1024.0) * Math.Sin(4.0 * latitudeRad)
                    - (35.0 * eSquared3 / 3072.0) * Math.Sin(6.0 * latitudeRad));

                double UTMEasting = UTM_F0 * n * (
                            A + (1.0 - t + c) * Math.Pow(A, 3.0) / 6.0
                            + (5.0 - 18.0 * t + t * t + 72.0 * c - 58.0 * ePrimeSquared) * Math.Pow(A, 5.0) / 120)
                        + 500000.0;

                double UTMNorthing = UTM_F0 * (M + n * Math.Tan(latitudeRad) * (A * A / 2.0
                                + (5.0 - t + (9.0 * c) + (4.0 * c * c)) * Math.Pow(A, 4.0) / 24.0
                                + (61.0 - (58.0 * t) + (t * t) + (600.0 * c) - (330.0 * ePrimeSquared)) * Math.Pow(A, 6.0) / 720.0));

                // Adjust for the southern hemisphere
                if (latitude < 0)
                    UTMNorthing += 10000000.0;

                return new UTM(UTMEasting, UTMNorthing, UTMZone, longitudeZone);
            }

            public void WGS84ToOSGB36() {
                var wgs84 = new RefEll(6378137.000, 6356752.3141);
                double a = wgs84.maj;
                double b = wgs84.min;
                double eSquared = wgs84.ecc;
                double phi = Deg2rad(this.Latitude);
                double lambda = Deg2rad(this.Longitude);
                double v = a / (Math.Sqrt(1 - eSquared * SinSquared(phi)));
                double H = 0; // height
                double x = (v + H) * Math.Cos(phi) * Math.Cos(lambda);
                double y = (v + H) * Math.Cos(phi) * Math.Sin(lambda);
                double z = ((1 - eSquared) * v + H) * Math.Sin(phi);

                double tx = -446.448;
                double ty = 124.157;
                double tz = -542.060;
                double s = 0.0000204894;
                double rx = Deg2rad(-0.00004172222);
                double ry = Deg2rad(-0.00006861111);
                double rz = Deg2rad(-0.00023391666);

                double xB = tx + (x * (1 + s)) + (-rx * y) + (ry * z);
                double yB = ty + (rz * x) + (y * (1 + s)) + (-rx * z);
                double zB = tz + (-ry * x) + (rx * y) + (z * (1 + s));

                var airy1830 = new RefEll(6377563.396, 6356256.909);
                a = airy1830.maj;
                b = airy1830.min;
                eSquared = airy1830.ecc;

                double lambdaB = Rad2deg(Math.Atan(yB / xB));
                double p = Math.Sqrt((xB * xB) + (yB * yB));
                double phiN = Math.Atan(zB / (p * (1.0 - eSquared)));
                for (var i = 1; i < 10; i++) {
                    v = a / (Math.Sqrt(1 - eSquared * SinSquared(phiN)));
                    double phiN1 = Math.Atan((zB + (eSquared * v * Math.Sin(phiN))) / p);
                    phiN = phiN1;
                }

                double phiB = Rad2deg(phiN);

                this.Latitude = phiB;
                this.Longitude = lambdaB;
            }

            public void OSGB36ToWGS84() {
                var airy1830 = new RefEll(6377563.396, 6356256.909);
                double a = airy1830.maj;
                double b = airy1830.min;
                double eSquared = airy1830.ecc;
                double phi = Deg2rad(this.Latitude);
                double lambda = Deg2rad(this.Longitude);
                double v = a / (Math.Sqrt(1 - eSquared * SinSquared(phi)));
                double H = 0.0; // height
                double x = (v + H) * Math.Cos(phi) * Math.Cos(lambda);
                double y = (v + H) * Math.Cos(phi) * Math.Sin(lambda);
                double z = ((1 - eSquared) * v + H) * Math.Sin(phi);

                double tx = 446.448;
                double ty = -124.157;
                double tz = 542.060;
                double s = -0.0000204894;
                double rx = Deg2rad(0.00004172222);
                double ry = Deg2rad(0.00006861111);
                double rz = Deg2rad(0.00023391666);

                double xB = tx + (x * (1.0 + s)) + (-rx * y) + (ry * z);
                double yB = ty + (rz * x) + (y * (1 + s)) + (-rx * z);
                double zB = tz + (-ry * x) + (rx * y) + (z * (1.0 + s));

                var wgs84 = new RefEll(6378137.000, 6356752.3141);
                a = wgs84.maj;
                b = wgs84.min;
                eSquared = wgs84.ecc;

                double lambdaB = Rad2deg(Math.Atan(yB / xB));
                double p = Math.Sqrt((xB * xB) + (yB * yB));
                double phiN = Math.Atan(zB / (p * (1 - eSquared)));
                for (var i = 1; i < 10; i++) {
                    v = a / (Math.Sqrt(1 - eSquared * SinSquared(phiN)));
                    double phiN1 = Math.Atan((zB + (eSquared * v * Math.Sin(phiN))) / p);
                    phiN = phiN1;
                }

                double phiB = Rad2deg(phiN);

                this.Latitude = phiB;
                this.Longitude = lambdaB;
            }

            public override string ToString() {
                return "'" + Math.Round(this.Latitude, 6).ToString(CultureInfo.InvariantCulture) + " "
                    + Math.Round(this.Longitude, 6).ToString(CultureInfo.InvariantCulture);
            }

            public DMS ToDMS() {
                //Input= xd(long) and yd(lat)
                //Output = xdd xm xs (long) and ydd ym ys (lat)
                double yd = Math.Abs(Latitude);
                double ydd = Math.Floor(yd);
                double ym = Math.Floor(60d * (yd - ydd));
                double ys = 3600d * (yd - ydd - ym / 60d);
                if (yd < 0)
                    ydd = -ydd;
                double xd = Math.Abs(Longitude);
                double xdd = Math.Floor(xd);
                double xm = Math.Floor(60d * (xd - xdd));
                double xs = 3600d * (xd - xdd - xm / 60d);
                if (xd < 0)
                    xdd = -xdd;
                return new DMS {
                    LongDegree = xdd, LongMinute = xm, LongSecond = xs, LongHemis = Longitude > 0 ? 'W' : 'E',
                    LatDegree = ydd, LatMinute = ym, LatSecond = ys, LatHemis = Latitude > 0 ? 'N' : 'S',
                };
            }
        }

        public class DMS {
            public double LatDegree;
            public double LatMinute;
            public double LatSecond;
            public char LatHemis;
            public double LongDegree;
            public double LongMinute;
            public double LongSecond;
            public char LongHemis;

            public static bool TryParse(object text, out DMS result) {
                result = null;
                if (!(text is string) || string.IsNullOrEmpty((string)text)) return false;
                var match = Regex.Match(((string)text).Trim(), "^([^°]+)° ?([^']+')? ?([^\"]+\")? ?([NS]) +([^°]+)° ?([^']+')? ?([^\"]+\")? ?([WE])$");
                double latD, latM, latS, lonD, lonM, lonS;
                if (!(match.Groups[0].Success
                    && double.TryParse(match.Groups[1].Value, NumberStyles.Float, CultureInfo.InvariantCulture, out latD)
                    && double.TryParse("0" + match.Groups[2].Value.TrimEnd('\''), NumberStyles.Float, CultureInfo.InvariantCulture, out latM)
                    && double.TryParse("0" + match.Groups[3].Value.TrimEnd('"'), NumberStyles.Float, CultureInfo.InvariantCulture, out latS)
                    && double.TryParse(match.Groups[5].Value, NumberStyles.Float, CultureInfo.InvariantCulture, out lonD)
                    && double.TryParse("0" + match.Groups[6].Value.TrimEnd('\''), NumberStyles.Float, CultureInfo.InvariantCulture, out lonM)
                    && double.TryParse("0" + match.Groups[7].Value.TrimEnd('"'), NumberStyles.Float, CultureInfo.InvariantCulture, out lonS)))
                    return false;
                var latCar = match.Groups[4].Value[0];
                var lonCar = match.Groups[8].Value[0];
                result = new DMS {
                    LatDegree = latD, LatMinute = latM, LatSecond = latS, LatHemis = latCar,
                    LongDegree = lonD, LongMinute = lonM, LongSecond = lonS, LongHemis = lonCar
                };
                return true;
            }

            public LL ToLL() {
                double latitude = Math.Abs(LatDegree) + LatMinute / 60d + LatSecond / 3600d;
                double longitude = Math.Abs(LongDegree) + LongMinute / 60d + LongSecond / 3600d;
                return new LL(LatHemis == 'S' ? -latitude : latitude, LongHemis == 'E' ? -longitude : longitude);
            }

            public string ToStringDMS() {
                return "'" + LatDegree.ToString(CultureInfo.InvariantCulture) + "°"
                    + LatMinute.ToString(CultureInfo.InvariantCulture) + "'"
                    + Math.Round(LatSecond, 2).ToString(CultureInfo.InvariantCulture) + "\"" + LatHemis
                    + " " + LongDegree.ToString(CultureInfo.InvariantCulture) + "°"
                    + LongMinute.ToString(CultureInfo.InvariantCulture) + "'"
                    + Math.Round(LongSecond, 2).ToString(CultureInfo.InvariantCulture) + "\"" + LongHemis;
            }

            public string ToStringDM() {
                return "'" + LatDegree.ToString(CultureInfo.InvariantCulture) + "°"
                    + Math.Round(LatMinute + LatSecond / 60, 6).ToString(CultureInfo.InvariantCulture) + "'" + LatHemis
                    + " " + LongDegree.ToString(CultureInfo.InvariantCulture) + "°"
                    + Math.Round(LongMinute + LongSecond / 60, 6).ToString(CultureInfo.InvariantCulture) + "'" + LongHemis;
            }

            public string ToStringD() {
                return "'" + Math.Round(LatDegree + LatMinute / 60 + LatSecond / 3600, 6).ToString(CultureInfo.InvariantCulture) + "°" + LatHemis
                    + " " + Math.Round(LongDegree + LongMinute / 60 + LongSecond / 3600, 6).ToString(CultureInfo.InvariantCulture) + "°";
            }
        }

        public class OSRef {
            public double Easting { get; private set; }
            public double Northing { get; private set; }

            // References given with OSRef are accurate to 1m.
            public OSRef(double easting, double northing) {
                this.Easting = easting;
                this.Northing = northing;
            }

            public LL ToLatLng() {
                var airy1830 = new RefEll(6377563.396, 6356256.909);
                double OSGB_F0 = 0.9996012717;
                double N0 = -100000.0;
                double E0 = 400000.0;
                double phi0 = Deg2rad(49.0);
                double lambda0 = Deg2rad(-2.0);
                double a = airy1830.maj;
                double b = airy1830.min;
                double eSquared = airy1830.ecc;
                double phi = 0.0;
                double lambda = 0.0;
                double E = this.Easting;
                double N = this.Northing;
                double n = (a - b) / (a + b);
                double M = 0.0;
                double phiPrime = ((N - N0) / (a * OSGB_F0)) + phi0;
                do {
                    M =
                      (b * OSGB_F0)
                        * (((1 + n + ((5.0 / 4.0) * n * n) + ((5.0 / 4.0) * n * n * n))
                          * (phiPrime - phi0))
                          - (((3 * n) + (3 * n * n) + ((21.0 / 8.0) * n * n * n))
                            * Math.Sin(phiPrime - phi0)
                            * Math.Cos(phiPrime + phi0))
                          + ((((15.0 / 8.0) * n * n) + ((15.0 / 8.0) * n * n * n))
                            * Math.Sin(2.0 * (phiPrime - phi0))
                            * Math.Cos(2.0 * (phiPrime + phi0)))
                          - (((35.0 / 24.0) * n * n * n)
                            * Math.Sin(3.0 * (phiPrime - phi0))
                            * Math.Cos(3.0 * (phiPrime + phi0))));
                    phiPrime += (N - N0 - M) / (a * OSGB_F0);
                } while ((N - N0 - M) >= 0.001);
                double v = a * OSGB_F0 * Math.Pow(1.0 - eSquared * SinSquared(phiPrime), -0.5);
                double rho =
                  a
                    * OSGB_F0
                    * (1.0 - eSquared)
                    * Math.Pow(1.0 - eSquared * SinSquared(phiPrime), -1.5);
                double etaSquared = (v / rho) - 1.0;
                double VII = Math.Tan(phiPrime) / (2.0 * rho * v);
                double VIII =
                  (Math.Tan(phiPrime) / (24.0 * rho * Math.Pow(v, 3.0)))
                    * (5.0
                      + (3.0 * TanSquared(phiPrime))
                      + etaSquared
                      - (9.0 * TanSquared(phiPrime) * etaSquared));
                double IX =
                  (Math.Tan(phiPrime) / (720.0 * rho * Math.Pow(v, 5.0)))
                    * (61.0
                      + (90.0 * TanSquared(phiPrime))
                      + (45.0 * TanSquared(phiPrime) * TanSquared(phiPrime)));
                double X = Sec(phiPrime) / v;
                double XI =
                  (Sec(phiPrime) / (6.0 * v * v * v))
                    * ((v / rho) + (2 * TanSquared(phiPrime)));
                double XII =
                  (Sec(phiPrime) / (120.0 * Math.Pow(v, 5.0)))
                    * (5.0
                      + (28.0 * TanSquared(phiPrime))
                      + (24.0 * TanSquared(phiPrime) * TanSquared(phiPrime)));
                double XIIA =
                  (Sec(phiPrime) / (5040.0 * Math.Pow(v, 7.0)))
                    * (61.0
                      + (662.0 * TanSquared(phiPrime))
                      + (1320.0 * TanSquared(phiPrime) * TanSquared(phiPrime))
                      + (720.0
                        * TanSquared(phiPrime)
                        * TanSquared(phiPrime)
                        * TanSquared(phiPrime)));
                phi =
                  phiPrime
                    - (VII * Math.Pow(E - E0, 2.0))
                    + (VIII * Math.Pow(E - E0, 4.0))
                    - (IX * Math.Pow(E - E0, 6.0));
                lambda =
                  lambda0
                    + (X * (E - E0))
                    - (XI * Math.Pow(E - E0, 3.0))
                    + (XII * Math.Pow(E - E0, 5.0))
                    - (XIIA * Math.Pow(E - E0, 7.0));

                return new LL(Rad2deg(phi), Rad2deg(lambda));
            }

            public override string ToString() {
                return "(" + this.Easting + ", " + this.Northing + ")";
            }

            public string ToSixFigureString() {
                var hundredkmE = Math.Floor(this.Easting / 100000);
                var hundredkmN = Math.Floor(this.Northing / 100000);
                var firstLetter = string.Empty;
                if (hundredkmN < 5) {
                    if (hundredkmE < 5) {
                        firstLetter = "S";
                    } else {
                        firstLetter = "T";
                    }
                } else if (hundredkmN < 10) {
                    if (hundredkmE < 5) {
                        firstLetter = "N";
                    } else {
                        firstLetter = "O";
                    }
                } else {
                    firstLetter = "H";
                }

                var secondLetter = "";
                var index = (int)(65.0 + ((4.0 - (hundredkmN % 5.0)) * 5.0) + (hundredkmE % 5.0));
                var ti = index;
                if (index >= 73) index++;
                secondLetter = ((char)index).ToString();

                var e = (int)Math.Floor((this.Easting - (100000 * hundredkmE)) / 100);
                var n = (int)Math.Floor((this.Northing - (100000 * hundredkmN)) / 100);
                var es = e.ToString();
                if (e < 100) es = "0" + es.ToString();
                if (e < 10) es = "0" + es.ToString();
                var ns = n.ToString();
                if (n < 100) ns = "0" + ns.ToString();
                if (n < 10) ns = "0" + ns.ToString();

                return firstLetter + secondLetter + es + ns;
            }
        }

        public class UTM {
            public double Easting { get; private set; }
            public double Northing { get; private set; }
            public int LongitudeZone { get; private set; }
            public char LatitudeZone { get; private set; }

            public static bool TryParse(string text, out UTM result) {
                result = null;
                if (string.IsNullOrEmpty(text)) return false;
                var split = text.Trim().Split(' ');
                var match = Regex.Match(split[0], @"^(\d+)([A-Z])$");
                double easting, northing;
                int zone;
                if (!(split.Length == 3 && match.Groups[0].Success
                            && int.TryParse(match.Groups[1].Value, NumberStyles.Integer, CultureInfo.InvariantCulture, out zone)
                            && double.TryParse(split[1].Trim(), NumberStyles.Float, CultureInfo.InvariantCulture, out easting)
                            && double.TryParse(split[2].Trim(), NumberStyles.Float, CultureInfo.InvariantCulture, out northing)))
                    return false;
                char emisphere = match.Groups[2].Value[0];
                try {
                    result = new UTM(easting, northing, emisphere, zone);
                    return true;
                } catch { }
                return false;
            }

            public UTM(double easting, double northing, char latZone, int lngZone) {
                this.Easting = easting;
                this.Northing = northing;
                this.LatitudeZone = latZone;
                this.LongitudeZone = lngZone;
            }

            public LL ToLL() {
                var wgs84 = new RefEll(6378137, 6356752.314);
                double UTM_F0 = 0.9996;
                double a = wgs84.maj;
                double eSquared = wgs84.ecc;
                double ePrimeSquared = eSquared / (1.0 - eSquared);
                double e1 = (1 - Math.Sqrt(1 - eSquared)) / (1 + Math.Sqrt(1 - eSquared));
                double x = this.Easting - 500000.0; ;
                double y = this.Northing;
                double zoneNumber = this.LongitudeZone;
                char zoneLetter = this.LatitudeZone;

                double longitudeOrigin = (zoneNumber - 1.0) * 6.0 - 180.0 + 3.0;
                //TO Check
                // Correct y for southern hemisphere
                if (zoneLetter < 'N')
                    y -= 10000000.0;

                double m = y / UTM_F0;
                double mu =
                  m
                    / (a
                      * (1.0
                        - eSquared / 4.0
                        - 3.0 * eSquared * eSquared / 64.0
                        - 5.0
                          * Math.Pow(eSquared, 3.0)
                          / 256.0));

                double phi1Rad = mu
                    + (3.0 * e1 / 2.0 - 27.0 * Math.Pow(e1, 3.0) / 32.0) * Math.Sin(2.0 * mu)
                    + (21.0 * e1 * e1 / 16.0 - 55.0 * Math.Pow(e1, 4.0) / 32.0)
                      * Math.Sin(4.0 * mu)
                    + (151.0 * Math.Pow(e1, 3.0) / 96.0) * Math.Sin(6.0 * mu);

                double n = a / Math.Sqrt(1.0 - eSquared * Math.Sin(phi1Rad) * Math.Sin(phi1Rad));
                double t = Math.Tan(phi1Rad) * Math.Tan(phi1Rad);
                double c = ePrimeSquared * Math.Cos(phi1Rad) * Math.Cos(phi1Rad);
                double r = a * (1.0 - eSquared) / Math.Pow(1.0 - eSquared * Math.Sin(phi1Rad) * Math.Sin(phi1Rad), 1.5);
                double d = x / (n * UTM_F0);
                double latitude = (
                  phi1Rad - (n * Math.Tan(phi1Rad) / r) * (d * d / 2.0
                    - (5.0 + (3.0 * t) + (10.0 * c) - (4.0 * c * c) - (9.0 * ePrimeSquared)) * Math.Pow(d, 4.0) / 24.0
                    + (61.0 + (90.0 * t) + (298.0 * c) + (45.0 * t * t) - (252.0 * ePrimeSquared) - (3.0 * c * c)) * Math.Pow(d, 6.0)
                        / 720.0)) * (180.0 / Math.PI);

                double longitude = longitudeOrigin + ((d
                    - (1.0 + 2.0 * t + c) * Math.Pow(d, 3.0) / 6.0
                    + (5.0 - (2.0 * c) + (28.0 * t) - (3.0 * c * c) + (8.0 * ePrimeSquared) + (24.0 * t * t)) * Math.Pow(d, 5.0) / 120.0)
                        / Math.Cos(phi1Rad)) * (180.0 / Math.PI);

                return new LL(latitude, longitude);
            }

            public override string ToString() {
                return "'" + this.LongitudeZone.ToString(CultureInfo.InvariantCulture) + this.LatitudeZone.ToString(CultureInfo.InvariantCulture)
                    + " " + Math.Round(this.Easting).ToString(CultureInfo.InvariantCulture)
                    + " " + Math.Round(this.Northing).ToString(CultureInfo.InvariantCulture);
            }
        }

        private class RefEll {
            public double maj;
            public double min;
            public double ecc;

            public RefEll(double maj, double min) {
                this.maj = maj;
                this.min = min;
                this.ecc = ((maj * maj) - (min * min)) / (maj * maj);
            }
        }

        static double SinSquared(double x) {
            return Math.Sin(x) * Math.Sin(x);
        }

        static double CosSquared(double x) {
            return Math.Cos(x) * Math.Cos(x);
        }

        static double TanSquared(double x) {
            return Math.Tan(x) * Math.Tan(x);
        }

        static double Sec(double x) {
            return 1.0 / Math.Cos(x);
        }

        static double Deg2rad(double x) {
            return x * (Math.PI / 180);
        }

        static double Rad2deg(double x) {
            return x * (180 / Math.PI);
        }

        /**
         * Take a string formatted as a six-figure OS grid reference (e.g.
         * "TG514131") and return a reference to an OSRef object that represents
         * that grid reference. The first character must be H, N, S, O or T.
         * The second character can be any uppercase character from A through Z
         * excluding I.
         *
         * @param ref
         * @return
         * @since 1.1
         */
        OSRef GetOSRefFromSixFigureReference(string value) {
            char char1 = value[0];
            char char2 = value[1];
            var east = int.Parse(value.Substring(2, 3)) * 100;
            var north = int.Parse(value.Substring(5, 3)) * 100;
            if (char1 == 'H')
                north += 1000000;
            else if (char1 == 'N')
                north += 500000;
            else if (char1 == 'O') {
                north += 500000;
                east += 500000;
            } else if (char1 == 'T')
                east += 500000;

            int char2ord = (int)char2;
            if (char2ord > 73) char2ord--; // Adjust for no I
            double nx = ((char2ord - 65.0) % 5.0) * 100000.0;
            double ny = (4 - Math.Floor((char2ord - 65.0) / 5.0)) * 100000.0;
            return new OSRef(east + nx, north + ny);
        }

        /**
         *  Work out the UTM latitude zone from the latitude
         *
         * @param latitude
         * @return
         * @since 0.2
         */
        static char GetUTMLatitudeZoneLetter(double latitude) {
            if ((84 >= latitude) && (latitude >= 72)) return 'X';
            else if ((72 > latitude) && (latitude >= 64)) return 'W';
            else if ((64 > latitude) && (latitude >= 56)) return 'V';
            else if ((56 > latitude) && (latitude >= 48)) return 'U';
            else if ((48 > latitude) && (latitude >= 40)) return 'T';
            else if ((40 > latitude) && (latitude >= 32)) return 'S';
            else if ((32 > latitude) && (latitude >= 24)) return 'R';
            else if ((24 > latitude) && (latitude >= 16)) return 'Q';
            else if ((16 > latitude) && (latitude >= 8)) return 'P';
            else if ((8 > latitude) && (latitude >= 0)) return 'N';
            else if ((0 > latitude) && (latitude >= -8)) return 'M';
            else if ((-8 > latitude) && (latitude >= -16)) return 'L';
            else if ((-16 > latitude) && (latitude >= -24)) return 'K';
            else if ((-24 > latitude) && (latitude >= -32)) return 'J';
            else if ((-32 > latitude) && (latitude >= -40)) return 'H';
            else if ((-40 > latitude) && (latitude >= -48)) return 'G';
            else if ((-48 > latitude) && (latitude >= -56)) return 'F';
            else if ((-56 > latitude) && (latitude >= -64)) return 'E';
            else if ((-64 > latitude) && (latitude >= -72)) return 'D';
            else if ((-72 > latitude) && (latitude >= -80)) return 'C';
            else return 'Z';
        }

    }


}
